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Abstract. We present the first stable dynamical numerical evolutions of the 
Einstein equations in terms of a conformally rescaled metric on hyperboloidal 
hypersurfaces extending to future null infinity. Axisymmetry is imposed in 
order to reduce the computational cost. The formulation is based on an earlier 
axisymmetric evolution scheme, adapted to time slices of constant mean curvature. 
Ideas from a previous study by Moncrief and the author are applied in order to 
regularize the formally singular evolution equations at future null infinity. Long- 
term stable and convergent evolutions of Schwarzschild spacctime are obtained, 
including a gravitational perturbation. The Bondi news function is evaluated at 
future null infinity. 



PACS numbers: 04.20.Cv, 04.20.Ha, 04.25.D-, 04.30.-w 

1. Introduction 

A major problem in numerical relativity that is currently attracting considerable 
interest is the global treatment of entire asymptotically flat spacetimes. This is 
relevant both to the modelling of the gravitational radiation emitted by compact 
astrophysical sources, and to several questions in mathematical relativity such as 
the nonlinear stability of black holes and the quantitative asymptotic behaviour of 
perturbations of such spacetimes. 

In the standard approach based on the Cauchy formulation, one truncates the 
spatial domain at a finite distance from the source, where boundary conditions must 
be supplied (see [T] for a review). Apart from being compatible with the Einstein 
constraint equations and yielding a well-posed initial-boundary value problem, these 
boundary conditions must carefully take into account the propagation of gravitational 
radiation. For instance, crude choices of boundary conditions are known to lead 
to incorrect results on radiation tails in the evolution of perturbed black hole 
spacetimes [2j |3| . Recently, improved absorbing boundary conditions that minimize 
spurious reflections of gravitational radiation have been derived and implemented 
[4j El [6] . While constituting a significant improvement, this approach still relies on the 
assumption that the Einstein equations may be linearized about a given background 
spacetimc so that incoming and outgoing radiation can be identified — an assumption 
that may not always be valid [7]. In the full nonlinear theory there is no satisfactory 



An axisymmetric evolution code for the Einstein equations on hyperboloidal slices 2 

quasi-local definition of a gravitational energy flux at a finite distance from the source 
[5] . This fact also complicates the extraction of gravitational radiation in an invariant 
manner. 

Gravitational radiation is only unambiguously defined at future null infinity 

[9] [10]. Thus it is very desirable to include ^ + in the computational domain. 
There arc several ways of doing this. One approach is characteristic evolution, whereby 
spacetimc is foliated by null hypersurfaces that can be compactified towards J'^ (see 
[TT] for a review). This works well in the far field but characteristic foliations are ill 
behaved in strong-field regions due to the formation of caustics in the null congruences 
generating the hypersurfaces. A compromise is Cauchy-characteristic matching, where 
a Cauchy evolution in the interior is matched to a characteristic evolution in an 
exterior region. In |12j this has recently been applied to the a posteriori extraction 
of gravitational radiation from a Cauchy evolution (with finite outer boundary) of a 
binary black hole coalescence. 

A more efficient solution is a foliation of spacetime into hyperboloidal 
hypersurfaces, which are everywhere spacelike but approach future null infinity rather 
than spacclikc infinity. Such surfaces are not Cauchy surfaces — the foliation covers 
only the region to the future of some initial hyperboloidal hypersurface but not, for 
instance, the region around spacelike infinity. This is not a serious disadvantage 
for the applications we have in mind as we are mainly interested in computing the 
gravitational radiation at ^ + and studying the late-time behaviour of perturbations. 
Hyperboloidal foliations have the advantage of being as flexible as standard Cauchy 
foliations in the interior, so that a variety of gauge conditions can be used that have 
been known to be successful in strong-field situations. In [13] Friedrich reformulated 
the Einstein equations for a conformally rescaled metric in such a way that they 
are manifestly regular at J? + . This symmetric hyperbolic system is larger than the 
actual Einstein equations, e.g. it contains evolution equations for the Wcyl curvature. 
For reviews of the theoretical development as well as numerical implementations, we 
refer the reader to |14[ [15] [16]. Perhaps most relevant to the present work is an 
axisymmetric reduction of the regular conformal field equations by Frauendiener and 
Hein |T7] that was applied to the numerical evolution of Minkowski spacetime and the 
boost-axisymmetric solutions of Bicak and Schmidt [18] . 

Despite these successes, most recent work in numerical relativity has used 
formulations of the Einstein equations (usually on truncated Cauchy hypersurfaces) 
rather than the regular conformal field equations. For this reason, it is worth 
investigating whether some formulation of the Einstein equations themselves can be 
adapted to hyperboloidal foliations reaching out to ,f + . As for the regular conformal 
field equations, one applies a conformal transformation to the metric in order to map 

to a finite coordinate location in a compactified coordinate system, where the 
conformal factor vanishes. When written in terms of the conformal metric, the Einstein 
equations contain inverse powers of the conformal factor, terms that are singular at 

. Numerical implementations must face the question of how to deal with these 
formally singular terms, which may be potential sources of instabilities. 

One proposal for solving the Einstein equations on hyperboloidal slices was 
developed by Zenginoglu |19| . It combines a conformal transformation of the spacetime 
metric (where the conformal factor can essentially be freely specified) with a certain 
choice of gauge source functions for the Einstein equations in (generalized) harmonic 
gauge [20) . Preliminary numerical results on spherically symmetric evolutions of 
Schwarzschild spacetime were reported in section 2.5 of [21]. However, it is not yet 
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clear how the formally singular terms at ^ + are to be evaluated numerically in more 
general situations. 

Here we follow a different approach due to Moncricf that is based on an ADM- 
like [35] formulation of the Einstein equations with a constant-mean-curvature (CMC) 
slicing condition. In [23] we showed explicitly how the formally singular terms in the 
evolution equations at can be evaluated there in terms of conformally regular 
geometric data. The constraint equations that occur in this approach were solved 
numerically in [24] and initial data on CMC slices describing various configurations of 
single and binary black holes were constructed. The present paper is concerned with 
the evolution problem. We assume spacetime to be axisymmetric in order to speed up 
the code so that different ways of treating J? + numerically can be experimented with. 
We apply the ideas of |23] to an earlier axisymmetric evolution scheme on Cauchy 
slices [25] (see also (26] [27] for related schemes). The objective of the paper is to show 
that long-term stable evolutions are possible, including gravitational radiation. 

The paper is organized as follows. In section [5] we explain our notation and 
gauge choices and obtain the constraint and evolution equations. Particular emphasis 
is placed on how the formally singular evolution equations can be regularized at 

. In section [3] we describe our numerical method for solving these equations. 
The numerical results are presented in section [4] As a first test problem we evolve 
Schwarzschild spacetime and demonstrate the convergence of the numerical solution. 
Next, a gravitational perturbation is included. A regular expression at ^ + for the 
Bondi news function [9] describing the outgoing gravitational radiation is derived and 
evaluated numerically. We conclude and discuss possible directions for future work in 
section [5] 

2. Formulation 

2.1. Definitions and gauge choices 

We assume that spacetime is axisymmetric with Killing vector £. For simplicity we also 
assume here that £ is hypersurface orthogonal. Spherical polar coordinates t, r, 8, <j> 
are chosen such that £ = d/d<p. 

The spacetime metric ^g a p is written as 

Wg a p = 4>- m 'Y a p, (1) 
where ip is the conformal factor and ' 4 '7 Q( g the conformal metric. Greek indices a, /3, . . . 
are spacetime indices (t,r,9,<fi). The conformal factor ip is positive and approaches 
zero at future null infinity, which we put at a fixed coordinate location r = 1. The 
conformal metric ^"/ a /3 is assumed to be regular there. 
We decompose the metric in ADM form as 

(4) 5Q ^dx Q da; /3 = ~a 2 dt 2 + g ij (dx i + p i dt)(dx j + ftdt), (2) 

where gij is the spatial metric, a the lapse function and f3 l the shift vector. Latin 
indices . . . are spatial (r, 9, <j)). Similarly, we write the conformal spacetime metric 
as 

[A) la pdx a dx fi = -& 2 dt 2 + ^j(dx l + /3 i dt)(dx j + (3 3 dt), (3) 
where %j = tp 2 gij is the conformal spatial metric and a = i/joe the conformal lapse. 
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The spatial coordinates are chosen such that ffi = and the conformal spatial 
metric takes the form 

- , ,(l.r'(b-' = e 2 " sine (dr 2 + r 2 d9 2 ) + r 2 sin 2 6»d0 2 . (4) 

This is sometimes called quasi-isotropic gauge. We note that this choice of gauge 
differs from the spatial harmonic gauge used in [23J ||| 
The extrinsic curvature is defined by 

d t gij = -2aKij + Cpgij, (5) 

where C denotes the Lie derivative. Preservation of the conditions ggg = r 2 g rr and 
g r e = under the evolution equation ([5]) implies 

(3 r , r - P\e - r- x p r = a(K r r - K\) = aU, (6) 

P r .e + r 2 p\ r = 2aK r e = 2&K r e . (7) 

Similarly as in [25j we have defined the following components of the extrinsic curvature, 

K r r -K\ = ^U, K r e = r - 2 K e r = ipk r g, (8) 

and, motivated by regularity concerns on the axis of symmetry [25] HH1 HZ] , we also 
set 

o — tj> = ip sin 9 W. (9) 

It is easy to see that the quantities U, K r g and W defined in ([5]) and © have the 
same conformal weight as the traceless momentum 

i^A'/i (10) 

used in [23], where g denotes the determinant of g^ . Note also that because K r g, W 
and U are either off-diagonal or differences of diagonal components of the extrinsic 
curvature, they can in fact be regarded as the three independent components of the 
traceless part of the extrinsic curvature in our case. 

As in |23] we require the time slices to have constant mean curvature, 

//'•*' A',, = -K = const. (11) 

Note the slightly awkward sign in the definition of the constant K; with our sign 
convention for the extrinsic curvature ([5j we require K to be positive so that the slices 
reach future null infinity. Preservation of (jllj) under the Einstein evolution equation 
for the extrinsic curvature yields the following elliptic equation for the conformal lapse, 

= a. rr + 2r~ 1 a r + r~ 2 (<5 00 + f <5 0) 

-^sa(ri, rT + r^rj.r + r- 2 r} M + 2r^ 2 ^r\ fi - r~ 2 ?y) 

- \^- 2 ae 2 ^K 2 



-fa 



P r (P r - 2 AO + r- 2 P e {Pe - 2A e 



_| ae 2 ^ \(U +\sW) 2 + \{sW) 2 +r- 2 {K r e ) 2 . (12) 
Here and in the following we use the abbreviations 

s = sin 8, c = cos ( 
upper-case Latin indices A, B, . . . ranging over r and 



Pb = ip 1 ip,B, A-b = a 1 a^ B , (13) 



\ It is, in a sense, its two-dimensional analogue: if we define a two-dimensional metric Hab by the 
part of the line element J2} orthogonal to the Killing vector, HAgdx A dx B = e 2r > BlnB (dr 2 + r 2 d9 2 ), 

then V c = H AB (T c ab — r° Ab) = 0, where Y c ab are the Christoffcl symbols of Hab an d Y° ab 
are the Christoffcl symbols of the flat metric {Hab with r\ = 0). [The author thanks V. Moncricf for 
pointing this out.] 



An axisymmetric evolution code for the Einstein equations on hyperboloidal slices 5 

2.2. Constraint and evolution equations 

As usual in ADM-like reductions, the Einstein equations split into constraint and 
evolution equations. The Hamiltonian constraint is 

= W + 2r"V,r + r- 2 (^. ee + §V,e) - [Wv) 2 + ^ 2 (<M 2 ] 

-\si\)(j] irr + r _1 ?/, r + r~ 2 r].gg + 2r~ 2 f n fi - r~ 2 i]) 

-\^ 2sr > h(U + \sW) 2 + \{sW) 2 + r- 2 {K r e ) 2 

+ i^-i e ^K 2 . (14) 

The momentum constraints are 

= C r = \U, r + \sW, r + r- 2 K r g fi + r" 2 ^(f + 2«7 + 2sn fi - 2Pg) 

+U{sr], r ~ f P- + 2T- 1 ) + sWir- 1 - §P r ), (15) 

= C e = -\U fi + \sW, e + K r e, r + 2K r g(sr], r - P r + r^ 1 ) 

+U{-c7 1 - sr?, e + \P B ) + W(±c - |sP fl ). (16) 
We have the following evolution equations, 

i), t = /3 r ^ r + - (f P° + r- 1 ^) - \&{K + ^{U + 2sW), (17) 

r),t = P r V,r + f3 e v,e + -J e V + (s^b - &W, (18) 

W t = (3 r W r + (3 e W e + (2§ f3 e + r^p^W 

+c -^r~ 2 [-( s - l a ie )j + 2aV" 1 (s-V,e),e] 
-ae- 2sri {n trr + 2r~ 1 r], r + r- 2 \r) fi g - rj + c^r])^ 

+(A r - 2P r )7 lr - r^s^iAg - 2P 9 )(sn, e + cr])} 
-2s- Y K r gfi\ r + ffis- 1 [sW^sVF + \U- i/j^K) + 3r- 2 (K r g ) 2 ] (19) 

K r 0>t = p r K r g, r + lfK r g fi + ZJ e k r e + e- 2s "(-a re + 2^"V,re) 

+fic" 2s " [(A r - 2P r + r-^iarifi + en) + (Ag - 2Pg){sji, r + r" 1 ) 

+a lr ] + l&k r e (sW - ip^K + 2U) - r 2 p\ r U, (20) 
U,t = P r U, r + P 6 Ufi + U{- s P e + r- l p r ) 

+c-' 2s,, [-a, rr + r- 2 a fie + 2a^~ l (tp , rr - r _2 ^,ee)] 
+5e" 2s, '[(i r - 2P r )(2s?7,r + r" 1 ) + 2sr~\ T 

-2r- 2 (i e -2P e )( S 7 ? , e + c7 7 +f)] 
+ iaf/(2siy - 2tP~ 1 K +U) + AK r e {f3\ r - ar- 2 K r g). (21) 

We observe that all the equations are manifestly regular at the axis of symmetry 
6 = when the parities of the fields with respect to 6 = are taken into account 
(table [T]). These can be inferred from the general behaviour of regular axisymmetric 
tensor fields [2"5] . 

2.3. Partially constrained evolution and hyperbolicity 

We determine the conformal factor ijj by solving the momentum constraint ()14[) 
rather than the evolution equation (|17[) . On the other hand, the extrinsic curvature 
variables W, K r g and U are evolved using the evolution equations (fT9|) -([2~T j) and the 
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Table 1. Angular parities of the fundamental variables: even (+) or odd (— ). 
The parities about 8 = follow from regularity on axis, those about 9 = n/2 from 
the reflection symmetry we impose. 



Variable 


a 


(3 r 


P 9 






w 




u 


6 = 


+ 


+ 




+ 








+ 


6 = tt/2 


+ 


+ 




+ 


+ 


+ 




+ 



momentum constraints (|15[) and (|16l) are not solved explicitly. Hence we adopt a 
partially constrained evolution scheme. 

Of course the constraints are preserved by the evolution equations, which imply 
the following subsidiary system, 

C\ t ~ i3 r C\ r + (3 e C\ e - |c«KCV + r- 2 C%), (22) 

C\ t ~ (3 r C\ r + f3 e C r .g + |e- 2s "aV _1 C*, r , (23) 

C e . r ~ (3 r C e . r + (3 C e ,g + |e- 2 "'c^- 1 C t )fl , (24) 

(25) 

where the symbol ~ signifies that only the principal parts have been displayed, and the 
remaining terms are homogeneous in the constraints. If the Hamiltonian constraint 
is enforced exactly, C* = 0, then the remaining evolution equations for C r and C are 
clearly hyperbolic (they are advection equations along the shift). If the Hamiltonian 
constraint is not enforced then the full constraint evolution system is not hyperbolic 
(the characteristic speeds 

A r = p r ± J^e-^ai, \ e = /?" ± J\ r-V'ai (26) 

are complex) . This might be cured by adding multiples of the Hamiltonian constraint 
to the evolution equations; however, we have not found such a constraint addition that 
both yields a hyperbolic constraint evolution system and is compatible with regularity 
of the main evolution equations on the axis of symmetry. Hence a free evolution 
scheme that solves none of the constraint equations during the evolution docs not 
appear to be feasible. However, our formulation includes three elliptic equations for 
the gauge variables a and (3 A anyway and solving one additional elliptic equation for 
the conformal factor is not overly expensive. 

If we assume that a, (3 A and tp are given by the solution of the elliptic equations, 
then the pair of evolution equations (|18|) and (|19jl forms a wave equation to principal 
parts, 

(fit - p A d A fi 1 ~ a 2 e- 2s "(77, rr + r-\gg), (27) 

and the evolution equations (|20p and ([21]) are advection equations along the shift to 
principal parts. Hence the evolution equations are hyperbolic. A rigorous proof of 
the well posedness of the present mixed elliptic-hyperbolic hyperboloidal initial value 
problem remains to be carried out. For a similar mixed elliptic- hyperbolic formulation 
on a Cauchy foliation, such a proof was given in in |29j . 

2-4- Regularity at future null infinity 

The CMC slicing condition (|12p . the constraint equations p^|) - ([TB| and the evolution 
equations for the extrinsic curvature (fT9|) - (|2~lT) are formally singular at ^ + , where 
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the conformal factor ip vanishes. Our elliptic solver (section 13.4)) is well adapted to 
such degenerate elliptic equations. However, the right-hand sides of the evolution 
equations cannot be evaluated at J? + in their present form. In [53] two methods for 
obtaining regular forms of the evolution equations at were presented: the first 
based on Taylor expansions, the second on an argument due to Penrose [30]. For the 
numerical results presented in section |4]wc followed the second approach but we have 
also implemented the first approach, as discussed further below. 

Penrose showed that provided certain smoothness assumptions hold (in particular, 
the conformal metric must be C 3 up to the boundary), the conformal Weyl tensor (i.e., 
corresponding to the conformal metric ^jafi) vanishes at . This implies [23] that 
the electric part of the physical Weyl tensor (i.e., corresponding to ^g a p) also vanishes 
at J+ , 

E al = n^n 5 C aM5 [^g] = 0, (28) 

where n a is the unit normal to the t = const slices and = denotes equality at . 
It is straightforward to compute the Weyl tensor for the metric ^'jap defined in ([3J 
and ([4]). Any time derivatives are substituted using the evolution equations. (Here a 
computer algebra program is very useful; we used REDUCE [31].) Setting E ai = 
immediately gives regular expressions for the formally singular terms in (fl"9]) (|2T|) . 
Thus we obtain the following manifestly regular evolution equations at ^ + , 

W t = (3 r W r + l3 e W, e + (2f /3 e + r^p^W - e^r" 2 ^- 1 ^ e ),e 

+ae~ 2s71 {q^ rr + 2r _1 ?7 ir + r~ 2 [r] M - Tj + c(s _1 77),e] 

-A r r]. r + r~ 2 s~ 1 A e (sn^ + crj)} 
-2s- l K r e p e ,r + 5[4 s -V- 2 (A^ e ) 2 - UW], (29) 

K r e , t = p r K r e ,r + f3 e K r e,e + £ s P e K r e - e~ 2s ^a, re 

+ae~ 2sr i [{A r - r- 1 )(s? / , e + CT)) + A g (sn, r + r" 1 ) - cr), r ] 
+2aK r e (sW + U)-r 2 l3 tr U, (30) 

U,t = P r U,r + P B U,e + U{^ e + r- l p r ) + e" 2s "(-a „■ + r~ 2 a M ) 

+& C - 2s,1 [A r (2 SVtr + r- 1 ) - 2r- 2 (Ag - § )(s Vfi + cq) - 2sr"V] 
+aU(2sW + U)+ AK r (p\ r - ar- 2 K r e ). (31) 

The other approach to regularity at discussed in [53] is based on an expansion 
of all the fields in (finite) Taylor series with respect to radius r about JP + . Substituting 
these expansions in the singular constraint equations, one obtains conditions on the 
fields and their radial derivatives that can be used in order to evaluate the formally 
singular terms in the evolution equations. This approach is more general in that it 
does not require Penrose's smoothness assumption. 

As necessary conditions for a regular evolution at , we obtained in [53] a 
set of conditions that had been discovered in [32] , namely that the shear of and 
the radial components 7r tr " of the traceless momentum vanish. In the axisymmetric 
formulation we consider here, these conditions read 

W = - c- s,1 r], r , U = ic- ST 'sry, r , K r g = 0. (32) 

One might be tempted to impose these conditions as Dirichlet conditions on the 
extrinsic curvature at ^ + , thus avoiding to evaluate the singular evolution equations 
there. However, numerical evolutions with these Dirichlet boundary conditions were 
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found to be unstable. This is perhaps not surprising because at J 2 "" 1 " all the fields are 
purely outgoing and hence imposing boundary conditions there is ill posed (although of 
course the conditions (j32")l are satisfied for a given solution of the field equations) . It is 
worthwhile to observe though, as we already did in |23j . that the evolution equations 
(|19[) - (j2"Tj) contain remarkable damping terms that cause potential violations of the 
regularity conditions (|32|) to decay exponentially in the neighbourhood of J? + : 

u t = l&i^Ku (33) 

for each u £ {W,K r g, U}. The term has the "right sign" because K > 0. 

Instead of imposing the regularity conditions (f3"2"j) directly, however, we can 
extract more information from the Taylor expansions that allows us to evaluate 
the singular terms in the evolution equations directly |23j . The resulting modified 
evolution equations at J !+ are detailed for the current axisymmetric scheme in 
|Appendix B| Their main difference to (|29| -(|3T |) lies in the presence of terms involving 
first radial derivatives of the extrinsic curvature at (cf. equation (|B.9jl ). We found 
this version of the evolution equations at J? + to be numerically stable as well, and 
the results agree with those in section [4] within numerical error (sec |Appcndix B] for 
details). 

It is worth noting here that the quasi-isotropic spatial gauge we use differs from 
the spatially harmonic gauge of [23]- Nevertheless, we were able to apply the same 
procedure for evaluating the formally singular terms at J? + . This provides some 
substance to the claim in [23] that the precise form of the spatial gauge conditions is 
inessential for the method to be applicable. 



3. Numerical method 



Having obtained a set of gauge conditions, constraints and evolution equations, we now 
describe how we solve them numerically. First we summarize the evolution scheme, 
i.e. which variables are evolved using which equations. Next we specify our choice for 
the computational domain and the discretization of the partial differential equations. 
We then focus on the two main aspects of the code, the time integration scheme and 
the elliptic solver. The code has been written in C++. Unless otherwise indicated all 
routines have been developed from scratch. 



3.1. Evolution scheme 

Our fundamental variables are a, (3 A ,ifj,n,W,K r g and U. The variables rj, W, K r e 
and U are evolved using the evolution equations (|18 p — (|2i p . At each time step we solve 
the CMC slicing condition (|T2"j) for a, the quasi-isotropic gauge conditions (JB]) and © 
for (3 A , and the Hamiltonian constraint (HU for ip. Details on how the quasi-isotropic 
gauge conditions are treated are given in |Appcndix A.l The momentum constraints 
are only solved initially ( |Appcndix A.2[ ). 



3.2. Domain and discretization 

For the applications considered in this paper, the spatial numerical domain is taken 
to be the quarter of an annulus, r m j n ^ r $J 1 and ^ # ^ tt/2. The inner boundary 
at r = r m j n will be taken to lie within a black hole event horizon (black hole excision) , 
and the outer boundary at r — 1 corresponds to . The axis of symmetry is at 
9 = 0, and we impose an additional reflection symmetry about the plane 9 = tt/2. 



An axisymmetric evolution code for the Einstein equations on hyperboloidal slices 9 

The fields typically have the steepest gradients close to the black hole and hence 
it is advisable to use a non-uniform radial grid in order to provide more numerical 
resolution towards the inner boundary. We introduce a quadratic map 

r(x) = Qx 2 + (1 - r min - Q)x + r min , (34) 

where ^ Q < 1 is a constant. This map satisfies r(0) = r min and r(l) = 1. We use 
a uniform grid in x, i.e. we place grid points at xi = i/N r , ^ i ^ N r , where N r is 
the number of radial grid points. Note in particular that there are grid points right at 
the inner boundary and at . In the angular direction we use a uniform staggered 
grid: Oj = f (j — \)/Ng, 1 ^ j ^ Ng. Two layers of ghost points are added on either 
side, corresponding to j = —1,0 and j = Ng + i, Ng +2 . Values at these ghost points 
are set according to the angular parities of the various fields, which are listed in table 
[T] E.g. for a field u that is odd about 9 — we set Uio = —Ui t i and Ui_-\ = — u.;2, and 
for an even field we set — ui_\ and iti,-i = u^. Here Uij denotes the value of u at 
the grid point with indices (i, j). 

The spatial derivatives are discrctized using fourth-order accurate finite 
differences. One-sided differences are used in the radial direction at the boundaries; 
everywhere else centred differences are used (see |Appendix C| for details). 

3.3. Time integration 

A fourth-order Rungc-Kutta method is used in order to integrate the evolution 
equations forward in time. At each substep of this method, the elliptic equations are 
solved as described below in section l3~4l Since all the characteristics point towards the 
exterior of the domain both at the inner excision boundary and at ^ + , no boundary 
conditions are imposed there — the evolution equations are discretized using one-sided 
finite differences. At J? + the evolution equations ([T9 |) -(f2T ]) for the extrinsic curvature 
are replaced with their regularized versions (|29| -(f3T |) . Kreiss-Oliger dissipation [33] is 
added in order to ensure stability \ Appendix C[ ). No modifications to the dissipation 
operators or special smoothing operations at have been required. 

3.4- Elliptic solver and boundary conditions 

The elliptic equations are solved by means of a multigrid method [3U E5] using Full 
Approximation Storage (FAS). Because of the singular nature of the equations at 
.y + (r = 1), we have found it necessary to use line relaxation in the radial direction, 
i.e. all points on a line 9 = const are solved simultaneously. With our fourth-order 
finite-difference discretization this requires solving a penta-diagonal linear system of 
equations, for which a fast algorithm exist (we use the routines bandec and banbks 
in [36]). For the Hamiltonian constraint (|14j) an outer Newton- Raphson iteration is 
employed for each such line solve; all remaining elliptic equations are linear. The lines 
9 = const are then traversed and updated in ascending order j = 1, 2, . . . , Ng (Gauss- 
Seidel relaxation in the ^-direction). We have found this method to be very robust 
and effective, leading to a reduction of the residual of about an order of magnitude per 
W-cycle. No problems with convergence at the singular boundary were encountered. 

The elliptic equations require boundary conditions at the inner r = r m j n and 
outer r — 1 boundary (in addition to the parity boundary conditions at 9 — 0,7r/2), 
which we shall now specify. 
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At the outer boundary (J" + ) we take the angular shift to vanish, (3 = 0. The 
conformal factor vanishes there by definition, ip = 0. Preservation of this condition 
under the evolution equation (|17[) yields a boundary condition for the conformal lapse, 

a= -e SI? /3 r . (35) 

We remark that these boundary conditions for the lapse and shift ensure that the 
normal n a to the time slices is null at . The radial shift f3 r appearing in (f3"S"|) cannot 
be chosen freely but is determined by the quasi-isotropic gauge conditions. After each 
Gauss-Seidcl relaxation sweep we integrate ([7]) along the boundary to determine /3 r 
there. The integration constant can be fixed e.g. by requiring that the average of 
(3 r along the boundary agree with the value obtained from the exact Schwarzschild 
solution that we consider in section 2] 

At the inner excision boundary, we also take (3 e to vanish, and we freeze a to 
the value corresponding to the exact Schwarzschild solution. However, the conformal 
factor ip cannot be chosen freely. This is because ip obeys an evolution equation with 
outgoing characteristic speed at the excision boundary, equation (|17[) . In order to 
determine the correct value of ip there, we keep a copy of ip that is evolved using (|17[) . 
At the end of each timestep this copy of ip is reset to agree with the solution of the 
Hamiltonian constraint. This procedure was found to be crucial for ensuring that the 
momentum constraints are preserved during the numerical evolution. 



4. Numerical results 



In this section, we apply our code to the evolution of a (perturbed) Schwarzschild black 
hole. First we focus on the unperturbed Schwarzschild solution and demonstrate that 
long-term stable and convergent evolutions can be obtained. Next we include a small 
gravitational perturbation. The outgoing gravitational radiation as represented by the 
Bondi news function is computed at J? + . 



Schwarzschild spacetime 

First we describe how the exact solution is computed, which is needed both for the 
boundary conditions of some of the elliptic equations (section I3.4|) and in order to 
compare with the numerical solution during the evolution. The Schwarzschild solution 
in CMC slicing was derived in [37l [38] . Its line element is 

ds 2 = - (l - — ) dt 2 + \df 2 - — dt df + f 2 (d6 2 + sm 2 6 d^ 2 ) 7 (36) 

V r J f f 

where 

m^il-^+a 2 ) 1 ' 2 , «(*0 = : X-^' (37) 

and M (mass), K (mean curvature) and C are constants. 

As in |24j we introduce a new radial coordinate r (the coordinate used in our 
evolution scheme) by requiring the spatial metric to be conformally flat in the new 
coordinates. Hence r obeys 

% - -TTY ^ 

dr rj{r) 

This ordinary differential equation is integrated numerically to high precision using 
Mathcmatica. The boundary condition is that r — > 1 as f — > oo. 
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The metric variables defined in section 12.11 arc then given by 

i>= r -, & = —, /3 r = - — , /3 e =0, i] = 0. (39) 
r r r 

The extrinsic curvature variables are found to be 

3(7 

[/ = - — , W = 0, K r e=0, (40) 

and the vector V A introduced in order to solve the momentum constraints initially 
( |Appendix A.2 1 is 

V r = ^, V e =0. (41) 

For the numerical evolutions presented in this paper we choose M = 1, K = 1/2 
and C = 2 (although we have obtained stable evolutions across a wide range of 
parameter values). The inner boundary is placed at r mm = 0.05, just inside the 
horizon at f = 2 r = 0.0635. 

We run the simulation with two different numerical resolutions (N r , Ng) = (64, 8) 
and (128, 16). The radial grid is nonuniform with Q = 3/4 in (|34| in both cases. 

The time step At is limited by the Courant-Friedrichs-Lewy condition 

. ( Ar A9\ 

At<min — - lfTT , 42 

Vh4l \ V ±\J 

where Ar is the grid spacing and v± are the incoming and outgoing characteristic 
speeds of the evolution system (fl~7|) - (j2"Tj) in the r-direction (similarly in the 9- 
direction), and the minimum is taken over all grid points. For our choice of parameters 
the time step is constrained by the outgoing radial characteristic speed at , 

2K 

Y 

for the Schwarzschild solution. The time step is taken to be about 3 /4 of the maximum 
allowed value computed in this way: for the lower resolution we use At = 0.06. 

During the evolution, we compute the difference of the numerical solution with 
respect to the exact solution for each of the evolved variables and take its discrete L2 
norm. This is divided by the norm of the exact solution for the corresponding variable 
whenever the latter does not vanish. Of these results a vector 2-norm (square root of 
the sum of the squares) is taken. The resulting quantity represents a measure of the 
overall relative error of the numerical solution as a function of time and is plotted in 
figure [TJ The error grows only linearly with time. It decreases by a factor of about 
15 on average as the resolution is doubled, close to the expected convergence factor 
of 2 4 = 16 for a fourth-order accurate scheme. The main contribution to the error 
comes from the region close to the black hole horizon where the fields have the steepest 
gradients. The numerical solution remains smooth all the way up to during the 
entire evolution. 



v' + = -f3 r + e- s '"a = — (43) 



4-2. Perturbed Schwarzschild spacetime 



Next, we include a small gravitational perturbation in our evolutions. This is done 
by taking the variable 77, which vanishes for the exact Schwarzschild solution, to be a 
Gaussian initially, 

(r - r c ) 2 ~ 



i] = A sin 9 exp 



2a 2 



(44) 
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Figure 1. Relative error of the numerical solution with respect to the exact 
solution for an evolution of Schwarzschild spacctime, computed as described in the 
text. Two different numerical resolutions are shown, (N r ,Ng) = (64,8) (dashed 
line) and (128, 16) (solid line). 



The factor sind has been included for parity reasons (tabled]). We choose A = 10 -4 , 
r c = 0.5 and a = 0.05. We stress that these initial data are evolved with the full 
nonlinear Einstein equations. The linearized problem has been studied extensively; 
in particular, numerical evolutions of the Bardeen-Press equation |39j and the Regge- 
Wheeler-Zerilli equations [40] have been performed on hyperboloidal slices. 

In order to demonstrate the consistency and convergence of our code, we monitor 
the momentum constraints, which are not solved explicitly during the evolution. For 
the unperturbed Schwarzschild solution, each term in the ^-momentum constraint (|16[) 
vanishes separately. This is not true for the r- momentum constraint (|15|) ; it vanishes 
only because of a subtle cancellation of the terms 



for the Schwarzschild solution. Each separate term in this equation is large (in 
particular close to the inner boundary) and the violation of the r-momentum constraint 
in the perturbed evolutions is dominated by the finite-differencing error of these 
terms for the Schwarzschild background. Therefore we have divided the r-momentum 
constraint by the norm of a typical term in (|45p (roughly 10 2 ) in order to give a better 
idea of the relative size of the constraint violation. The results are shown in figure 
[2] The violation of the momentum constraints decreases by just over an order of 
magnitude as the resolution is doubled, reasonably close to the expected convergence 
factor of 2 4 = 16. 

Of prime interest is the gravitational radiation emitted by the perturbed black 
hole as measured at ^ + . It can be represented by the Bondi news function [5]. This 
function was defined originally by studying the asymptotic behaviour of the metric in 
a special coordinate system. This coordinate system will in general not agree with the 
coordinates used in a numerical evolution. Fortunately, Stewart [H] has developed 
a method for computing the news in an arbitrary coordinate system, working as we 
do here in terms of a conformally related metric. We begin by specifying a Ncwman- 
Pcnrose [12] null tetrad (I, n, m, fh) that is adapted to ,f + in the sense that 



-|P r +2r- 1 ) =0 



(45) 



m a d a i/> = 0. 



(46) 
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Figure 2. Momentum constraints (115P — (116P as functions of time for two different 
resolutions (N r ,Ng) = (64,8) (dashed lines) and (128,16) (solid lines). The r- 
momcntum constraint (left panel) has been normalized by the size of a typical 
term in (1451 ) as explained in the text. 



For definiteness we choose 

"-^(•-'l+'-O- (47) 

Apart from the requirement (|46|) the tetrad is arbitrary. The resulting expressions will 
be invariant under spin and boost transformations. From the conformal spacctime 
metric ' 4 '7 a/ s specified in ^ and (jH), we compute the conformal Ricci tensor ^R a p. 
The time derivatives in this expression are substituted by using the evolution equations 
(|17[) — (f2T|) . The Bondi news function is now given by 

N = #02 = m a fh 0W R a0 . (48) 

The resulting expression is not yet regular at but it can be written in a regular form 
by using the regularity conditions Q32p. the boundary condition (|35[) and expressions 
(|B.2[) and (|B.3[) for the radial derivatives of the conformal factor. Ultimately we obtain 

N = -e- 2sr > [s{r), rr + ri fi e) - s 2 (^ r ) 2 + cn,e - s -1 ?/] -e^cT 1 /3%, r - e - s W r .(49) 

The Bondi news function has spin weight —2 and hence we may expand it in spin- 
weighted spherical harmonics |43j _2Yfm with m = because of the axisymmetry, 

oo 

N(t, = N *(t) -2Y m {6). (50) 

1=2 

The I = 2 contribution is plotted in figure [3J The higher-^ contributions are 
considerably smaller (by at least a factor 10~ 4 ) and have not yet converged for the 
numerical resolutions used here. 

After an initial peak when the outgoing part of the wave reaches , we expect 
the news to consist of a ringing phase due to the quasi-normal mode excitations of 
the black hole, followed by a polynomial decay ("tail") due to the backscatter of 
the radiation off the curved background spacctime. These features are apparent in 
figure [3J However, whereas the ringing phase has converged quite well, the tail is still 
highly resolution-dependent. Considerably higher resolutions and/or a higher order 
of accuracy of the discretization will be needed in order to study the tail behaviour. 
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Figure 3. £ = 2 contribution to the Bondi news function l|49jl . evaluated at , 
as a function of time. Numerical results for two different resolutions are shown, 
(N r ,N e ) = (64, 8) (dashed line) and (128, 16) (solid line). 



We expect the quasi-normal mode ringing to have the form 

N e occ~ Ket sm(uJit + ip e ). (51) 

The decay rate K£ and frequency u>t are estimated by performing a \ 2 fit of (|51[) to 
the numerical data for the news at in the time interval 60 ^ t ^ 120. We find 
n, 2 = 0.0893 ± 0.0009 and w 2 = 0.3738 ± 0.0012. These values correspond to the 
high-resolution run and the errors have been estimated by comparing with the low- 
resolution run. This numerical result is consistent with the linearized-theory result 
from Leaver's continued-fraction method [44], K 2 = 0.08896 and UJ2 = 0.37367. 

5. Conclusions 

The objective of this paper was to demonstrate that stable numerical evolutions of the 
Einstein equations on hyperboloidal slices extending to future null infinity ^ + can be 
obtained. While there has been extensive numerical work |14 [ I15 [ [T6 ] using the regular 
conformal field equations |13j on such foliations, this is, as far as the author is aware, 
the first result using directly the Einstein equations in a dynamical situation. The 
advantage of this approach is that it is more similar to formulations of the Einstein 
equations on (truncated) Cauchy hypcrsurfaccs that have been used successfully in 
numerical relativity for a long time. 

In this paper, we assumed spacetime to be axisymmetric in order to be able 
to experiment with different ways of treating numerically in a relatively short 
computational time. The axisymmetric reduction of the Einstein equations we used 
(section[2]) is similar to previous works [25J [26l [27] but with a constant-mean-curvaturc 
slicing condition instead of maximal slicing. 

The main new difficulty in adapting a formulation of the Einstein equations 
on Cauchy slices to hyperboloidal slices is the treatment of the formally singular 
terms in the evolution equations at .y + . We applied two different methods [23] 
in order to obtain regular versions of these equations at (section |2.4[) . For the 
evolutions presented in this paper, we used Penrose's |30j result on the vanishing of 
the conformal Weyl tensor at J !+ . For this to hold, certain smoothness assumptions 
on the conformal metric at must be satisfied. An alternate method that does 
not require these assumptions is based on Taylor-expanding the fields near and 
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evaluating the singular elliptic equations order by order; with this information, the 
apparently singular terms in the evolution equations can be evaluated ( |Appcndix B[ ) . 
This alternate form of the evolution equations at was found to be stable as well 
and to yield comparable results to within numerical error. Future work should focus 
on this latter approach as it requires less restrictive smoothness assumptions and is 
most likely compatible with the so-called polylogarithmic behaviour at J" + [32] . 

Our numerical implementation (section[3]) uses fairly standard techniques: fourth- 
order finite differences, the method of lines with a fourth-order Rungc-Kutta scheme 
for the evolution equations, and a multigrid elliptic solver. The only non-standard 
modification to the multigrid solver is the use of a line relaxation in the radial direction 
in order to deal with the singular elliptic equations. With this method, no problems 
with the convergence of multigrid at the singular boundary were encountered. For the 
evolution equations, the usual Kreiss-Oliger dissipation {33j sufficed to ensure stability; 
no special smoothing operations at J? + were required. 

We were able to evolve Schwarzschild spacetime for 1000M (M being the black 
hole mass) without any signs of instabilities and with nearly perfect fourth-order 
convergence (section |4|. Next we included a small gravitational perturbation. In 
order to demonstrate the consistency and convergence of the scheme we monitored 
the momentum constraints during the evolution. We derived a regular expression 
for the Bondi news function at J? + HI] representing the outgoing gravitational 
radiation. From this we computed the frequency and decay rate of the quasi-normal 
mode radiation and found good agreement with the tabulated values from linearized 
theory. The numerical resolutions we were able to reach were not yet high enough in 
order to study the subsequent polynomial decay (the "tail"). 

Currently the code is too slow in order to make runs with higher resolutions 
feasible (the run at the higher resolution in section 14.11 took about two weeks on a 
single processor). Not much effort has been spent on optimizing the computational 
efficiency of the code and there is certainly some room for improvement here. It 
seems likely, however, that a finite-difference method of higher than fourth order or 
a pseudospectral method will need to be used in order to reach the computational 
accuracy needed to study tails. 

Although our multigrid solver has optimal computational complexity (linear in the 
number of unknowns) [341 135] , the fact that our formulation contains elliptic equations 
is a major drawback on the computational speed. This is likely to be even more severe 
in the case without spacetime symmetries. From this point of view, a hyperbolic 
formulation such as the one by Zcnginoglu [T!5] seems advantageous. It would be 
interesting to see whether the methods of [23] can be applied to that formulation as 
well so that regular versions of the evolution equations at are obtained. 

In this first study we imposed axisymmetry in order to be able to experiment 
with different numerical treatments of ^ + more quickly. The next step will be the 
extension to spacetimes without any symmetries. As far as the numerical stability of 

is concerned, we do not expect any new difficulties here because unlike spherical 
symmetry, axisymmetry already shares with the non-symmetric case what we believe 
are the key challenges, namely the existence of a dimension tangential to the conformal 
boundary and the presence of gravitational radiation. It remains to be seen which 
formulation of the Einstein equations will be more successful for this purpose. Wc 
expect the ideas and methods developed in this work to be applicable to a variety of 
formulations and gauge choices. 
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Appendix A. Further details on the solution of the elliptic equations 

Appendix A.l. Quasi-isotropic gauge conditions 

The quasi-isotropic gauge conditions ([6]) and ([7]) are formally equivalent to the 
(inhomogeneous) Cauchy-Ricmann equations. In order to solve them, one usually 
combines derivatives of these equations to obtain two decoupled Poisson equations, 
which can easily be solved. However, care must be taken so that the solution to the 
second-order Poisson equations is in fact also a solution to the original first-order 
Cauchy-Ricmann equations (in general the two may differ by integration constants). 
We use the following approach. 

The quasi-isotropic gauge conditions (JSJ) and (0 are 

= S r = f3 7 \ r ~ /3% - r^(5 r - aU, (A.l) 

= S e = p r , e + r 2 f3 e t r-2&K r g. (A.2) 

Derivatives of these are combined to obtain a scalar Poisson equation, 

= -S r ,g + r- 2 S e , r = p e , rr + r' 1 $\ r + r~ 2 f3 6 fi e - 2r- 1 &, r K r e + r- 2 a fi U 

- 2r~ 2 aK r g r + r~ 2 alJ fi + Ir^aK^. (A.3) 

This equation is solved for (3 e using multigrid, with homogeneous Dirichlet boundary 
conditions. 

However, we do not solve the corresponding Poisson equation for j3 r . Instead, we 
observe that once (3 e is known we can immediately integrate (|A.1[) and (|A.2|) to find 
f3 r . First, (|A.2|) is integrated along the outer boundary J? + to determine (3 r there. 
The integration constant can be chosen arbitrarily, e.g. by requiring the average of /3 r 
along J? + to agree with its value for the exact Schwarzschild solution. Next, (|A.1[) is 
integrated from each point on the outer boundary towards the interior. In this way we 
ensure that S r = everywhere. We have also enforced S e = at the boundary. What 
remains to be shown is that S e = in the interior. But this follows immediately 
from (|A.3|) : given that we already have S r = everywhere, this equation reduces 
to S 6 , r = 0, and together with the initial condition S e = this implies S e = 
everywhere. 

Apart from ensuring that the first-order conditions (|A.1[) and (|A.2[) are indeed 
satisfied, we also avoid solving another elliptic equation (for j3 r ) and thus save 
computational time. 

Appendix A.2. Momentum constraints 

For the intial data only, the momentum constraints are solved using a method similar 
to the one due to York [45] . We introduce a vector V A and write the extrinsic curvature 
variables as 

U = ip 2 {V\ r - V\g - r^V") = ?/> 2 VL, (A.4) 
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K r e = \^{V r s + r 2 V\ r ) = ^ 2 V + , (A.5) 

W = ip 2 W. (A.6) 

Note the similarity of (|A.4| and (|A.5|) with © and ©. The field W can be 
freely specified. The substitutions (|A.4j) - (|A.6[) completely decouple the momentum 
constraints (|15|) and (|16[) from the remaining elliptic equations: 

= \V r ,„ + \t- 2 V t m - \V\ re - \V b \e + \aW, r + r^sW 

+V-{sri >r + fr" 1 ) + r- 2 V + (s7ie + cr, + §§), (A.7) 
= \V\ r0 + \r 2 V\ rr + \V\ee - f r^V.e + ±5^,0 + f cVF 

-F_(sr/, e + 07) + 7+(«»;, r + 2r- : ). (A.8) 

We take F A to agree with the corresponding vector computed from the exact 
Schwarzschild solution at the radial boundaries, in particular we have V A = at J+ , 
and V e also vanishes at the inner boundary. The angular parity of V A is identical 
with that of f3 A (table [J). 

Appendix B. Alternate regularized evolution equations 

In this appendix we derive alternate regular forms of the evolution equations at 
.y + using the Taylor expansion method developed in [23j . Within a hyperboloidal 
hypcrsurface t — const, all the fields are expanded in finite Taylor series about 
(r = 1 in our case), 



u(x % ) = V —,uJ6, 4>)(r - 1)", u„ = lim d"u. (B.l) 
£ — ' n\ r/l 

n=0 

These series are substituted in the singular elliptic equations and evaluated order by 
order. Again we have found a computer algebra programme [3T] to be very helpful 
here. From the Hamiltonian constraint (fT4l) we obtain 



Vv = - \Kc sr i, (B.2) 
^, rr = -Ke"i(\sri, r + \), (B.3) 

^ rrr = Ke sv {ff(sr?),e + |(«/), 9 fl - § [(s?y),e] 2 + fs^rr + i|(^r) 2 + §s?7,r} ■ (B.4) 

The slicing condition (fT2"j) yields 

<5. r = <5(is77 jr + 1), (B.5) 
a, rr = a [ff (s77),e + \{sri).ee + |s7?,rr - \{sri^ r ) 2 + sn. r + l] 

^Ts^fi + l^fie ~ %&fi{srj)fi. (B.6) 
From the momentum constraints (| 15[) and (|16p we obtain, respectively, 

sW, r = -2f7 ir + |e- s "(s7 ?ir ) 2 > ( B - 7 ) 

iT^ = e~ sr ' [-CT/, r - |(«»/,r),fl + «7,r (*»/),*] ■ ( B - 8 ) 

Substituting the Taylor expansions (|B.lj) in the right-hand sides of the evolution 
equations (fT5)) - (f2"l"j) and using the above results (|B.2|) (|B.8|) for the radial derivatives 
of the fields, together with the regularity conditions (|32|) and the boundary condition 
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Figure Bl. £ = 2 contribution to the Bondi news function l|49jl at ^+ for an 
evolution using regularization at based on Taylor expansions (solid line) as 
compared with the approach based on Penrose's smoothness assumption (dashed 
line) that was used in section|4] The numerical resolution is (N r , Ng) = (128, 16). 



(|35[) . we obtain the following regular forms of the evolution equations at ^ + . 

W t = c~ 2sv {a [c^- 1 ^ + i] fiB + r], rr - s{r,_ r ) 2 - rj\ + a,e(f»7 + V,e) 

-{s-^ele} + e- S7 >aW, r + e^^Mcr? - 2f + sq fi ) - v, r e], ( B -9) 

K r e,t = 0, (B.10) 

U t = - \sW.t. (B.ll) 

Notice in particular how (|B.10[) and (|B. 1 1|) are compatible with the regularity 
conditions (|3"2"j) . 

Figure IB1I shows the Bondi news function at for a numerical evolution using 
the above Taylor-regularized evolution equations at J?^ . The result agrees well with 
the evolution based on the Penrose regularization technique that was discussed in 
section 12.41 and used for the evolution shown in figure |3l The discrepancy for times 
t > 150 can be attributed to numerical error — figure |3] shows that the solution is still 
highly resolution-dependent at those late times. The corresponding curves in figures 
[T] and [2] virtually overlap for the two regularization techniques. 

Appendix C. Finite-difference operators 

In order to compute derivatives with respect to r, we first compute finite-difference 
approximations to derivatives with respect to the map coordinate x and then apply 
the chain rule, 

_ dx _ f dx\ 2 d 2 x 

dr \ dr J dr* 

Here dx/dr and d 2 x/dr 2 are computed analytically from the given form of the mapping 
r(x), equation 

Fourth-order accurate finite difference operators are used with respect to both 
x and 9. It suffices to display their one-dimensional forms here. The centred finite 
difference operators are 

U 'i ~* 7Ih(- Ul - 2 ~ 8u i~l + 8u i+l - (C.2) 
U i ~~* VlK I (~ Ul -' 2 + 16m »-1 - 30u i + l 6u »+l _ u i+2), (C.3) 
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where h denotes the (uniform) grid spacing and the indices refer to the grid points 
where the quantities arc evaluated. The forward differences are 

u' -> ^(-251*0 + 48ui - 36u 2 + 16u 3 - 3it 4 ), (C.4) 

u' x — » i57r(^3iio — 10ui + 18u 2 — 6u 3 + u 4 ), (C.5) 

u o -> iw( 45u o _ 154 "i + 214u 2 - 156u 3 + 61u 4 - 10u 5 ), (C.6) 

u" -> 3^(10^0 — 15tti — 4u 2 + 14u 3 — 6U4 + U5), (C7) 

and expressions for the backward differences follow by symmetry. 

In the x-direction we use centred differences for 2 ^ i ^ N r — 2, forward differences 
at i = 0,1 and backward differences at i = N r , N r — 1. In the 0-dircction centred 
differences are used at all interior grid points 1 ^ j ^ Ng (note that ghost points are 
provided at j = -1, and j = N g + 1, Ng + 2). 

The expression (s _1 u) e is discretized by finite-differencing the function s~ 1- it 
with respect to 9. No such special differencing is used for (s u which we 
simply expand into s~ 1 u t gg — s~ 2 cu,g and then discretize in the standard way. Mixed 
derivatives u tX g are discretized by applying the finite-difference operators with respect 
to x and 9 subsequently in the obvious way. 

We use Kreiss-Oliger dissipation (33], whereby the operator 

(Qu)i = (^h 5 D 3 + D 3 _u). 

= ^(^-3 - 6m 4 _ 2 + 15u;_i - 20u.; + 15u.;+i - 6u. l+2 + u i+3 ) (C.8) 

is applied to each variable u both in the x and the 9 direction and the results added 
to the discretized right-hand side of the evolution equation for u. Here D± denote the 
first-order forward and backward finite-differencing operators, and a typical choice 
of the dissipation parameter is e = 1/2. In the x-direction the above dissipation 
operator is added at grid points 2 $J i ^ N r — 2 and in the ^-direction at grid points 
2 ^ j ^ Nq — 1. Including j = 2, Ng — 1 has been found to be essential in order to 
eliminate an angular instability. 
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